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Abstract 

In the present paper a lattice Boltzmann scheme is presented which ex¬ 
hibits an increased stability and accuracy with respect to standard single- 
or multi-relaxa-tion-time (MRT) approaches. The scheme is based on a 
single-relaxation-time model where a special regularization procedure is 
applied. This regularization is based on the fact that, for a-thermal flows, 
there exists a recursive way to express the velocity distribution function 
at any order (in the Hermite series sense) in terms of the density, velocity, 
and stress tensor. A linear stability analysis is conducted which shows 
enhanced dispersion/dissipation relations with respect to existing models. 
The model is then validated on two (one 2D and one 3D) moderately high 
Reynolds number simulations (Re ~ 1000) at moderate Mach numbers 
(Ma ~ 0.5). In both cases the results are compared with an MRT model 
and an enhanced accuracy and stability is shown by the present model. 


1 Introduction 


The lattice Boltzmann method (LBM) is a widely used tool for numerical sim¬ 
ulations of fluid flows. It has become over the years one the of the major engi¬ 
neering tools for computational fluid mechanics. It describes the flow thanks to 
the time evolution of the velocity distribution function which is only modified 
through the effect of inter-particle collisions. 

The most commonly used lattice Boltzmann collision model is the single re¬ 


laxation time model or BGK (for Bhathagar, Gross and Krook, see Bhatnagar 
et al. 1954] ). This model is able to asymptotically represent weakly compress¬ 
ible fluids (through a Ghapman-Enskog expansion, see Ghapman and Cowling 
[1960]). Nevertheless it suffers from stability issues especially at high Reynolds 


numbers. These issues are due to the “ghost-modes” (see Dellar|| 2001] for a 


discussion) which are non-physical moments present in any LBM simulation in 
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excess of the density (pressure), velocity and stress. This issue has been ad¬ 
dressed by several authors and several solutions have been proposed. The first 


is the multiple-relaxation-time (MRT) approach (see 

d’Humieres 

1992 , 

Lalle- 

mand and Luo 

20001, 

d’Humieres et al. |2002|, Dellar 2001 , 

Xu and Sagaut 


2011 , Xu et al. 2012 among others) which uses a more complex collision 


model involving several relaxation times adjusted with the help of a linear sta¬ 
bility analysis in order to optimize their dispersion/dissipation relations. The 
entropic approach ensures the positivity of the distribution functions and hence 
the unconditional stability by adding an i/-theorem to the BGK model. The 
effect of the i/-theorem is essentially to increase locally the viscous dissipation 


of the model (see Ansumali and Karlin 2002 , 

Boghosian et al. |2003 , 

Ghikata- 

maria et al. 2006 , 

Malaspinas et al. 

2008 among others). The regularization 

approach (see 

Latt and Ghopard 2006 , Zhang et al. 120061) which can be in- 


terpreted as a subclass of MRT methods where the “ghost-modes” are relaxed 
towards zero with characteristic time 1. Finally the selective viscosity models 
proposed by Ricot et al. [200^ use non-local low-pass filters to remove high 
frequency oscillations (which are responsible for the numerical instabilities) in 
order to increase the stability. 

In this paper we will first show a recursive way to compute the moments of 
the distribution function as long as the Chapman-Enskog expansion is valid (low 
Knudsen number) for the BGK collision operator. This recursive relation will 
then be used to “regularize” the distribution function and provide a very stable 
and accurate scheme even at moderately high Reynolds numbers and (relatively) 
high Mach numbers (smaller than one though). We also show that the present 
model is more accurate and more stable than the existing MRT methods by 
performing a linear stability analysis and several numerical benchmarks. 

The paper is structured as follows. In Sec.j^a reminder of fundamentals for 
fluid flows with the Boltzmann-BGK equation is presented. Then in Sec. the 
new model is proposed and analyzed. It is validated in Sec. [5on a 2D and a 3D 
benchmark. Finally the present work is concluded in Sec.^ and perspectives 
are given. 


2 The hydrodynamic limit of the BGK equation 


The following section aims at introducing the basic notations as well as showing 
the fundamentals of the expansion leading from the continuous Boltzmann-BGK 
equation to the Navier-Stokes equations. More details can be found in |Shan| 
et al. 2006 and Malaspinas [2009| for example. 


The Boltzmann equation describes the time evolution of the velocity density 
probability distribution of finding a particle with velocity ^ at position 

X and time t in terms of particle collisions only, and reads in absence of a force 
as 

+ (I • V)fix,tt) = m), (1) 


where D is the collision operator. Assuming also that the fluid is athermal 
(absence of temperature), the macroscopic fields of interest, the density p, the 
velocity u, and the stress tensor P are given by the following moments of the 
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distribution function 


p = ^ 


(2) 

pu = 


(3) 

p = ^ 

dc ccf{x,^,t), 

(4) 


where c = ^ — u is the microscopic velocity in the co-moving frame and cc 
denotes the tensor product of c with itself. 

The most widely used model for computational fluid dynamics for the col¬ 
lision operator is the BGK, single relaxation time approximation, in which the 
Boltzmann equation reads 

dtf{x,^,t) + (I • V)/(x,|,t) = , (5) 

where r the relaxation time, and is the local Maxwell-Boltzmann equilib¬ 
rium distribution function, which in non-dimensional units is given by 


fiO) 


(27r)^/2 


exp 




( 6 ) 


D being the physical dimension. 

Since we are interested in numerically solving Eq. ([^ in an efficient fashion 
that nevertheless represents accurately fluid flows, certain simplifications will be 
made. In particular instead of considering the complete form of the Maxwell- 
Boltzmann distribution function, only a polynomial approximation will be used. 

Following an idea by Shan et ^ 2006| (or Grad| |1949b| for the original use 
of this expansion in the frame of the Boltzmann equation), one can expand the 
distribution functions / and in Hermite polynomials up to an arbitrary 
order N (see [Grad 1949a for a summary on Hermite polynomials) 


N N 

77, 77, 


n—0 


n—0 


where the colon symbol stands for the full index contraction. The Hermite 
polynomials of order n and the associated Gaussian weight are noted and 
w{C) = exp(—^^/2) respectively. The Hermite coefficients of / and of de¬ 
gree n are respectively given by and a!'^\ From now on, we will always 
omit the superscript N and assume that the distribution function (and its equi¬ 
librium counterpart) is represented by its approximate form in terms of Hermite 
polynomials up to an arbitrary order N except when explicitly stated otherwise. 
The equilibrium coefficients can be easily computed and are found to be up to 
order three 


a 

a 


(0) 

0 

( 1 ) 


OCK 


a 


( 2 ) 

0a:/3 


a 


(3) 

Oa.^'y 


= P, 

(8) 

— ■) 

(9) 

— P'^a.^pi 

(10) 

- pUlQ/^U^Ury. 

(11) 
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In order to recover the macroscopic equations of motion related with the 
BGK equation, one must take moments of Eq. ([^. By taking the moments 
related with density (order zero) and momentum (order one) of this equation, 
one gets after some algebra and the use of Eqs. (H-Q and (l8])-@ 


dtp + V ■ {pu) = 0, 

dt{pu) + V • (puu) + V • P = 0. 


( 12 ) 

(13) 


These equations are obtained under the sole assumption of mass and momentum 
conservation (f = f = 0), or in other terms 


J (/-/^°^)=0, 
J d||(/-/(°))=0. 


Mass conservation (14) 

Momentum conservation (15) 


The momentum conservation equation still needs to be closed (a constitutive 
equation must be found for P). In order to do so, one can use the Chapman- 
Enskog expansion (see Chapman and Cowling [1960| , Huang 1987| ). Since the 
expansion in Hermite series is used for discretization purposes we will discuss 
the Chapman-Enskog expansion in this frame (although the Hermite series is 
not a prerequisite for performing the Chapman-Enskog expansion). 

The Chapman-Enskog expansion is based on the assumption that the dis¬ 
tribution function / is given by the sum of the equilibrium distribution, 
plus a small perturbation noted 




(16) 


where the equilibrium distribution is assumed to be given by Eq. 0. The 
perturbation, ~ 0{Kn) is of the order of the Knudsen number, 

Kn. As for / and one can express in terms of a Hermite series 




(17) 


n—0 


where is the Hermite coefficient of at order n. The derivation which 
is presented hereafter is not the standard one found in the literature and rather 
follows Huang [1987 


Replacing the Chapman-Enskog Ansatz in Eq. ([^, one obtains at the lowest 
order 

+ V)/(°) = --/W. (18) 

T 

Taking the zeroth and first order moments of this equation and using the mass 
and momentum conservation constrains on each equation respectively (f = 
= 0), one gets the inviscid Euler equations for mass, momentum and 
energy conservation 


dtp + V ■ (pu) = 0, (19) 

dt{pu) V • (puu) -\- Vp = 0, (20) 

where p = p is the perfect gas law (remember that there is no temperature). 
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The stress tensor can be decomposed in its Chapman-Enskog counterparts 


p^pio) ^p{i) =p7 + p(i)^ 


( 21 ) 


where = f ccf^^'> for j = 0,1 (j corresponding to the Chapman-Enskog 
index). Thus we are left with the computation of pl^^ which for simplicity is 
computed through the Hermite expansion of the distribution function. Let us 
dehne a\ ' the Hermite coefficient of order n of the off-equilibrium distribution 
function and express in terms of these Hermite coefficients 


p(l) _ „(2) 
.( 0 ) 


( 22 ) 


where we used that by construction a)'"'' = = 0. Then projecting Eq. ( |18[ ) 

on the Hermite basis, it follows that 


a a, 


t“o 


V • a 


(n-l-l) 


Va 


(n-l) 


perm I = - a\ , 


(23) 


where “perm” stands for all the cyclic index permutations. For n = 2 this 
equation becomes 


+ V • 


Vai 


( 1 ) 


, 1 ( 2 ) 

- perm I =- Ui , 


dt {puu) + V • {puuu) -I- {y{pu) + {y(pu))'^) =- 


1 (2) 


(24) 


By using Eqs. (|19[)-(20) to eliminate the time derivative terms, this equation 


2009 ) as 


can be rewritten (after some tedious algebra that can be found in Malaspinas 

^(2) ^ p(i) ^ -2tpS, 


where 


5=-(V.- 


[Vuf) . 


(25) 


(26) 


By comparing Eq. (25) with the Navier-Stokes equations, the transport coef¬ 
ficient p can be identified with the relaxation time through the following relation 

p = pr. (27) 

Finally the equations of motion obtained are the weakly compressible athermal 
Navier-Stokes equations 


dtp+V ■ (pu) = 0, 

dt{pu) + V • {puu) = — Vp -I- V • {2pS). 


(28) 

(29) 


3 Hierarchy of non-equilibrium moments and reg¬ 
ularization scheme 

In this section a novel theoretical approach is proposed for the athermal Boltzmann- 
BGK equation. A recursive formulation for non-equilibrium moments is shown 
to exist and a regularization technique is proposed for the discrete lattice Boltz¬ 
mann method. 
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3.1 Recursive properties of high order moments 

The particular structure of the moments of the equilibrium distribution in ab¬ 
sence of temperature allows for an elegant formulation of the high order (higher 
than two) non-equilibrium moments. 

The Hermite coefficients of order n of the equilibrium distribution can be 
recursively expressed as 


An) 


= ttg" ^'^u, and a® = p. 


(30) 


Using this relation and Eqs. (19)-(20l one can show that for n > 3 (see [A| for 
the proof) 




. ..OCn—l ‘ 




.( 2 ) 


'Oin-2^1,an- 


ia„ + Perm(a„)) , (31) 


where “perm(a„)” stands for all the cyclic index permutations of indexes from 
ai to an-i (an is never permuted). One therefore notices that a Hermite 
coefficient of order n can be expressed in terms of the the velocity and the 
Hermite coefficients of order order two and n — 1. This property allows to 
reconstruct the populations up to any order by only knowing its second order 
coefficient and the macroscopic velocity. 


3.2 Discretization of the microscopic velocity space 

We notice that only the Hermite coefficients of the distribution function are 
used in the Chapman-Enskog expansion. Therefore in order to asymptoti¬ 
cally recover the Navier-Stokes equations there is no need to use the complete 
Maxwell-Boltzmann equilibrium distribution but only a polynomial expansion 
of it. In order to discretize the velocity space one will use a Gauss-Hermite 
quadrature. The aim of this discretization is to exactly evaluate the integral of 
polynomials of order m with Gaussian weight as a sum 

='^w,pm(ii), (32) 

d i=0 


where and {^i}i=o 9 constant weights and abscissae 

respectively. 

In order to obtain asymptotically the weakly compressible limit of the BGK 
equation only polynomials of order m = 5 need to be integrated exactly. The 
associated most common quadratures (see Shan et al. 2006| ) for this case are 
given by the D2Q9 (in 2D) and the D3Q15, D3Q19, and D3Q27 lattice^(in 3D). 
These quadratures allow the definition of a set of velocity discretized distribution 
function noted as {fi}1Zo = {/(*j • In other terms, on each position 

X at time t one defines q independent values {/i}i=o- These quantities can 


therefore be represented on a g-dimensional basis (see work of d’Humieres 1992 


Dollar 2003 among others). We emphasize here is that there are two different 


spaces that must be distinguished: the velocity-discretized q-dimensional space 
and the d-dimensional physical space. 


'^The DdQg notation denotes a lattice of dimension d and with q quadrature points 
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The above quadratures only allow the exact representation of the populations 
up to second order in Hermite polynomials, one usually truncates Eq. Q to 
order two. This means that the equilibrium distribution is represented on a six¬ 
dimensional basis in 2D and respectively on a 10-dimensional basis in 3D, while 
it is living in a 9 (for the D2Q9 lattice) or 15, 19, or 27 (for the D3Q15, D3Q19, 
or D3Q27 lattices) dimensional space (depending on the quadrature used). 

The BGK equation discretized in microscopic velocity space then reads 

dth + ■ V)/, =-I f^°\p,u)) . (33) 


As pointed out above, the set populations {fi and /f°^) live in a g-dimensional 
space which means that they can be represented on a g-dimensional basis. In 

2^ 


d’Humieres 1992 , d’Humieres et al. 


or 


Dellar 2003 different bases are 


proposed for the expansion of the complete populations but never used to ex¬ 
press the equilibrium distribution. Here we propose to expand the equilibrium 
population on a complete basis as well. To this aim we will use an interesting 
property of the D2Q9 and the D3Q27 lattices, which is that the complete 9- and 
27-dimensional bases can be expressed in Hermite polynomials, a property that 
does not hold for the D3Q15 and D3Q19 quadratures. This property makes 
them particularly appealing in order to reuse all the calculations performed in 
the previous section. The distribution function is therefore written as (in 2D) 


=w, (^p + 

+ — a(3) -I-—^ 

y'^txxy^xxy ' ’ ^txyy^xyyj ' ' ^txxyy^xxyyJ * 

Respectively the equilibrium and off-equilibrium parts are expanded as 


f(0) (1 , 1 tj(2) 

f, =w,p ( 1 -b -b : uu 


2c® '^ixyy^x^yj “f ^^g'^ixxyy^x^y 


,(3) 


) 4c? 


/( 4 ) 


•I i y2c'^ ' 2c® y^ixxy^l,xxy ' ^ixyy^l,xyy 


—-wW U) 

4^8 ’^ixxyy^l,xxyy j ■ 


(34) 


(35) 


(36) 
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In 3D the equivalent expressions are given by 


fi =Wi 


ii ■ (pu) 


• a(2) 


— a(3) +77^^^ +77^^^ 

2^6 ixxy^xxy ' '^ixxz^xxz ' '^ixyy^xyy ' '^ixzz^xzz 

_l_ (,(3) j_ q(3) _|_ 277^^^ 

' '^iyzz^yzz ' '^iyyz^yyz ' ' ^ixyz^xyz j 


1 

4c® 


a(4) _|_ -^(4) ^(4) _|_ 

^txxyy xxyy ' ' ^za? 7 C 22 tctczz ' 


77^^) a(4) 

'^lyyzz yyzz 


+2 (77i!L.a('‘) 


I 77-"^^ + 77^"^^ 

'ixyzz xyzz ' ' ^ixyyz xyyz ' ' ^ixxyz xxyz 


— ^77^^^ 
4^10 y^ixx: 


(5) 


'ixxyzz^xxyzz 

77^®) a(6) 

g^l 2 ^txxyyzz^xxyyzz 


- 77^^^ 

' ^ixxyyz^xxyyz 


.(5) 


-T/(^) ^ 

^ixyyzz^xyyzz 


.(5) 


) 


(37) 


/f ^ =1"*^ ( 1 + 


^i-u 1 

2c4 


77 


( 2 ) . 


: ixtt 


^ i^fly'^luy + n^tlzUlu. + 'UfJyyU^ul + nfj,,u^ul 

~^'^\yzz'^y'^z + '^iyijz'^y'^z + ^^Lyz'^^a; Wy Mz^ 

— f77^‘‘^ ^2^2+77^^^ u2m2+77(4) ^2 2 

4^8 *xa:yy X ' '^ixxzz x^z ' ' ^lyyzz y z 

i^ixyzz'^xUyUz + '^ixyyz'^x'^y'^z + 'HixxyzUxUyUz^^ 

1 fn/{5) 2 2 I 0/(5) 2 2 I 0/(5) 2 2\ 

q^lO y~^ixxyzz^x^y^z “b ^ixxyyz^x^y^z. + ^ixyyzz^^^y^zJ 

(38) 


I ^0/(6) 2 2 2 

-r/• 7/ 7/ 7/ I 

g^l2 ^^2)2)7/1/22 ^7C*^2 i 5 


r(i) 

fl =Wi 


1 0 /( 2 ) . _( 2 ) 


— ( 
2c6 V 


'1/(3) ^(3) 

^ixxy^l^xxy 

i'^(3) (3) 

^iyzz^l.yzz 


-^(3) ^(3) 

^12)2)2^1,2)2)2 


_-^(3) ^(3) 

^ixyy^l,xyy ^ 12 ) 22 ^ 1 , 2)22 


_nji^) ^(3) -1-27-/^^^ 

^iyyz^l,yyz ' ^^ixyz^l,xyz 


^ -1-77 

4c® \^^^2)2)yi/^l,2)2)yy ^ 


(4) _(4) ^(4) (4) 

12)2)22^1,2)2)22 ^ ^iyyzz^l,yyzz 


+2 (77^!L,a^'‘^ 


4-77^^^ 4-77*'^^ 

12 ) 1 / 22 ^ 1 , 2)^22 ^ixyyz^l,xyyz ' ^ixxyz^l,xxyz 


1 

4cl3 




45 ) (5) 

"13) 2)y 22^1,2)3)^22 


^ixxyyz l,xxyyz 


^i3)yy22 l,2)i/y22 


.J_'W(6) „ 

g^l 2 ^ixxyyzz^l,xxyyzz 


,(6) 


( 39 ) 



The Hermite coefficients of the equilibrium distribution are the ones obtained 
from of the continuous equilibrium distribution for both the 2D and 3D cases, 
simplifying the computations. Another important remark is that not all the 
Hermite polynomials are used at each order. This is due to the fact that the 
quadrature is not accurate enough to represent exactly all the Hermite polyno¬ 
mials, but only the ones that are used in the formulas above. 


3.3 Chapman—Enskog expansion of the model 

The Chapman-Enskog expansion of this model asymptotically leads to the fol- 

( 2 ) 

lowing constitutive equation for ' in three dimensions 

„(2) 


( 2 ) 

( 2 ) 

A,yz 


2cIptSxx +Tdx{pul), 

(40) 

* 

2Cg pT Sxy) 

(41) 

P'^^xz ) 

(42) 

‘^clprSyy + Tdy{pul), 

(43) 

* 

PTSyz, 

(44) 

2cIptSzz +Tdz{pul), 

(45) 


where the terms are 0(Ma^) order error terms which are not present in 
the continuous case (see Eq. (p^). In the case of the order two standard BGK 
model (where is only expanded up to order two in Hermite polynomials) it 
reads 


^l,xx — ‘^CgPrSxx + '’'{dx^pv-x) ~^dyipUyUx) dz^pUzUx)), 


i,xy = -‘i.clprSxy +T{dx{puluy) -b dyipuxul) + dzipUxUyUz)), 




^1,XZ ‘^^sP'^^XZ “t“'r((9y “b Oxi^pUz^x) ^ z{,P^x‘^z)^ ^ 


— ‘^^sP'^'^yy “b '^i^yiP^y) ~^9x{pUxUy) -b dz{pUzUy)), 


Jl 2jCgPTSyz ~\~^{^dx^pUxUyUz') “b ^y^P^z^y) “b ^z^P'^y^z^^ 


«pL = -‘^ciprSzz + ridzipUx) +dx{puxui) + dy{puyui)), 


(46) 

(47) 

(48) 

(49) 

(50) 

(51) 


where the terms are the terms that are not present anymore in the new 
model. One can see that in our case the non-diagonal terms are exact. 

Furthermore since we now expand the distribution function up to a lim¬ 


ited order in Hermite polynomials the relations of Eq. (311 are not anymore 
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exactly verified. Nevertheless the difference between the exact relation and the 
error committed is compatible with the low compressibility approximation of 


the scheme (low Mach number approximation). In other terms Eq. (31) reads 
in the discrete case 




dn-l) 


+ Ua, ...Uo 


( 2 ) 


perm(a„)^ 


■0(Ma”+O, 


(52) 


where 0 stands for the order of the error committed. These terms are spurious 
terms that are due to quadrature (discretization in velocity space) errors in the 
expansion and are assumed to be small because they are one order higher in 
physical velocity. For a more detailed expression for the 0(Ma"~'’^) terms see 

m 


3.4 Time-space discretization 


The time space discretization of Eq. ( |33| ) is done as usual by integrating it 
along characteristics with the trapezoidal rule and making the following change 
of variables (see Debar |2001| for example) 


h = h^YT ~ ■ 

This leads to the following lattice Boltzmann method scheme 


(53) 


Mx + ^i,t + l) = f,{x,t)- - 

T 






(54) 


where r = t + 1/2. From now on the “bar” is omitted for brevity in the 
notations. 


3.5 Regularization scheme 


The numerical scheme used here is given by Eq. ( 

54 

) where fi is “regularized” 

at each iteration as done in 

Latt and Chopard 

2006 



(55) 


where is computed with Eqs. (36) or (39) depending on the dimension of 
the physical space. Furthermore the Hermite coefficients of Eqs. (36) and (39) 
(a^"\ with n > 2) are evaluated thanks to the recursive formulation of Eq. (31). 
For efficiency reasons Eq. (54) can thus be rewritten as 


f^{x + + 1) = 



(56) 


This is the novel scheme proposed in this paper. The model will be validated 
in Sec. 1^ and a comparison with an existing multiple-relaxation-time model will 
be performed. 
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3.6 Von Neumann linear stability analysis 


The aim of this section is to perform the linear stability analysis of the scheme. 
More details about the Von Neumann stability analysis of the lattice Boltzmann 
method can be found among others in Lallemand and Luo 2000 , Ricot et al. 


[2009| , Xu and Sagaut 2011 , Xu et al. 


2012 


By decomposing the distribution function into the sum of a stationary part 
(noted fi which must not be confused with the / of Subsec. 3.41 and a small 
fluctuating part, noted /', 

/. = /* + /', (57) 


1 


f(i) 


and by defining as the r.h.s. of Eq. (56) 

= /f ^ + (l - 

the linearized lattice Boltzmann scheme is found to be given by 
f'{x + ^„t+l) = ^Ay/j(a;,t), 


(58) 


(59) 


where A is defined as 


A _ ^ 

" dfj 


(60) 


fj=fj 


Working in Fourier space and looking for plane wave solutions Eq. (59) becomes 

//(fc,t+i) = ^A-%fc/afe,t), ( 61 ) 

j,k 


where = Sjk exp {—ik ■ $,j), with Sjk the Kronecker symbol and i = 

The eigenvalues, Aj, of the matrix allow to obtain the dispersion- 

dissipation relations of the numerical scheme, tOj{k) = zlogAj. While an an¬ 
alytic approach is used to determine these eigenvalues in [Lallemand and Luo| 
here we simply used a linear algebra package to determine numerically 


2000 


these eigenvalues. 

The Von Neumann stability analysis of the Navier-Stokes gives the following 
dispersion-dissipation relations 


Re(a;±) = ±||fc||cs-I-fc • M, (62) 

Im(a;±) = -||fc|p-(^—+ (63) 

Re(ws) = k ■ u, (64) 

=-WkW^v, (65) 


where rj = for the BGK model. 

We now compare the eigenvalues of the collision operator of the present 


model we the ones obtained for the MRT model proposed by Lallemand and 


Luo 2000 . The stability analysis depicted in this section has been performed 
for T = 0.5001 and u = (0.2,0) (corresponding to Ma = 0.346). As one can 
see from Fig. the dispersion relations for the present scheme are relatively 
similar to the ones obtained with an MRT approach except for the shear mode 
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Figure 1: Dispersion with respect to kx {ky = 0) for the present model and the 
MRT model with the Ux = 0.2, for the D2Q9 lattice, and r = 0.5001. This 
circled region highlights the region where there is a significant improvement of 
the dispersion relation of the present model as compared to the MRT model. 




Figure 2: Dissipation (right) with respect to kx {ky = 0) for the present model 
and the MRT model with the Ux = 0.2, for the D2Q9 lattice, and r = 0.5001. 
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Figure 3: Maximal value of the dissipation, max^ ujj{k), for the MRT (left) and 
present (right) models with Ux = 0.2, for the D2Q9 lattice, and r = 0.5001. 
The solid line represents the isocontour where max^ u!j{k) = 0. 


where the dispersion relation remains very close to the Navier-Stokes result 
until kx = TT which is not the case for the MRT (see the circled region of 
Fig. 0. For the dissipation (see Fig. 0 one clearly sees that for ky = 0 there 
is no unstable mode for the present model (Im(wj) < 0,\/kx,j) while this is not 
the case for the MRT model. Furthermore one can also notice that while the 
acoustic modes are only weakly dissipated, except for Wg, (and therefore the 
proposed scheme could be very suitable for aeroacoustic simulations) the other 
spurious modes are dissipated very fast and therefore the scheme should suffer 
of less numerical instabilities. The increased linear stability of the model is 
depicted in Fig. 0 where one can see that the imaginary part of the eigenvalues 
of the numerical scheme is always negative and therefore the scheme has an 
unconditional linear stability. This feature is not present in the case of the 
MRT model as for some values of k the imaginary part of the eigenvalues of the 
evolution operator become positive. Of course as one increases the magnitude 
of the velocity unstable modes will start to appear. The unstable modes start 
to appear at Ux = 0.248 (corresponding to Ma = 0.43) for t = 0.5001 as shown 
on Fig. 0 


3.7 Boundary conditions 


The aim of this subsection is not to give a detailed view of the way to implement 
boundary conditions since this topic is extensively treated in the literature (see 


Latt et al. 

2008 , 

Malaspinas et al. 2011|, Zou and He |1997 , 

Inamuro et al. 

1995 , 

Skordos 

1993 among others for some references). It rather explains how 


the proposed regularization model is compatible with all the cited boundary 
conditions. 

One way to deal with Dirichlet boundary conditions in the lattice Boltzmann 
method is to use the regularized procedure described in Latt et al. 2008 for 


example. The generic idea is to impose a velocity Ubc at the boundary. First one 
uses the symmetries of the lattice to compute p. Then it is possible to compute 
(which only depends on p and itbc)- Finally can be computed using 
a finite difference scheme through S or by using the symmetries of the lattice 
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Figure 4: Maximal value of the dissipation, maxj u!j{k), for the MRT (left) and 
present (right) models with Ux = 0.248, for the D2Q9 lattice, and r = 0.5001. 
The solid line represents the isocontour where max^ 0 Jj{k) = 0. 


(see 


Latt et al. 


2008 


). Then is computed by using the following formula 


f(i) _ . p(i) 


We notice that this formula is actually exact if Mbc = 0. By replacing Mbc by 
zero in Eq. (36) or (39) one is simply left with the equation above. Then if 
Mbc ^ 0, the procedure is exactly the same as for the two strategies discussed 
above and we should simply use Eqs. (36) or (39) for the computation of 


Finally one simply replaces all the populations on the boundary mesh point 


with the regularization formula (55). 


3.8 Differences with existing stabilization techniques 

Apart from the MRT-LBM there exists different techniques to increase the sta¬ 
bility and accuracy of the BGK-LBM scheme. In this subsection we discuss the 
major distinctions between the present scheme and some of these approaches. 
We will limit the discussion to the regularization, entropic, and selective viscos¬ 
ity filter techniques. 


3.8.1 The regularized model 


The existing class of regularized models (see Zhang et al. 2006 , Latt and 


Chopard 2006|) belongs to the same family as the model presented here. The 


general idea is the same as the one used for the present model. The collision 


operator is the same as Eq. (56) 


U{x + + 1) = + { 1 - ^ ] 


( 1 ) 


but instead of using the equilibrium distribution of Eqs. (35) (in 2D) or ( [3^ (in 
3D), and off-equilibrium distribution of Eqs. (36) (in 2D) or (39) (in 3D), one 


14 





































0.006 


0.005 
0.004 
^ 0.003 
S 0.002 
0.001 
0 


O 

o o Regularized 
X X Present 



- 0.001 



1.5 



2 2.5 3 


Figure 5: Dissipation (right) with respect to kx {ky = 0) for the present model 
and the MRT model with the Ux = 0.2, for the D2Q9 lattice, and r = 0.5001. 


truncates the series at order two in Hermite polynomials, which amounts to use 


=WiP ( 1 + 


■ u 


1 


/( 2 ) . 


+ -nr : nu 


n{l) _ — t qfy 

A 2d * 


w 


'i ^( 2 ) . (2) 


( 66 ) 

(67) 


This regularization technique removes all the moments of order higher than two 
in Hermite polynomials from the distribution function. These are considered 
as negligible in the asymptotic limit of the weakly compressible Navier-Stokes 
equation. The removal of these higher order terms affect the accuracy of the 
constitutive equation for the stress tensor as shown in Eqs. (40)- 


(45) and Eqs. (46)-(51). The present regularization not only provides a more 


accurate constitutive equation for the deviatoric stress but also preserves the 
recursive relation of the (for n > 3) terms. These differences lead to a major 
difference for the linear stability analysis of the two schemes. A comparison for 


of the dissipation (see Subsec. 3.6) for Ux = 0.2 and r = 0.5001 is shown in 
Fig. 0 One can see that while the uj± eigenvalues have very similar values for 
both models, the difference lies in the uig eigenvalue. For the regularized model 
this eigenvalue is positive (and therefore an unstable mode exists) while it is 
negative for the present model. Therefore one expects the present model to 
exhibit a much more stable behavior. 


3.8.2 The entropic model 


The entropic model (see among other 

Ansumali and Karlin 

2002 , iBoghosian 

et al. 

2003 , Chikatamarla et al. 

2006 

) is based on a different philosophy for 


the construction of the numerical scheme. The major difference is the existence 
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of an i?-function defined as 

^ = ( 68 ) 

i 

The assumption is then made that there exists a discrete _ff-theorem which 
states that 

1. The equilibrium distribution, minimizes the H function under the 

constrains that X)* = P and 

2. The H function is monotonically decreasing. 

This second condition is imposed through the following collision operator 

/(a; + |*,t+l) =/i-^ , (69) 

where a > 0 is computed such that 




(70) 


The collision operator of the entropic model guarantees an unconditional stabil¬ 
ity of the scheme. This comes nevertheless at a high computational cost since 
at each point and at each time the above non-linear implicit equation must 
be solved. As shown in several references (see Malaspinas et al. 2008 among 


others) the presence of the a parameter has as an effect to locally increase the 
viscosity (and therefore the dissipation). Therefore one would expect that the 
dissipation of the entropic scheme would be more important and therefore less 
suitable for acoustic propagation for example. 


3.8.3 The selective viscosity model 


In the selective viscosity model proposed by Ricot et al. |2009| the basic idea is 
to use the standard BGK-LBM collision operator (see Eq. (|54[)) and to increase 
the stability of the model by applying a low-pass hlter on the fi at each point 
and at each time step. The filtering operation is defined as 


D 


N 

dnMx + nej) 


(71) 


j = l n— — N 


where Bj are the D unit basis vectors of the Cartesian coordinate system, the dn 
are the 2N +1 filter coefficients, and a G [0,1] is the strength of the hlter. This 
hltering operation removes the high frequency oscillations responsible for the 
numerical instabilities appearing in the model. The different hlters proposed 
involve non-local computations (the hlters have width 


Ricot et al. 2009 


between three and nine mesh points) that impacts greatly the computational 
efficiency of the scheme since not only more operations must be performed at 
each mesh point but also the amount of communications (which are crucial for 
parallel efficiency) is also increased. In the model presented here no such non¬ 
local operations are performed reducing the computational cost with respect to 
the selective viscosity models. 
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Furthermore the filtering operation implies the existence of a cutoff which 


removes the high wavenumber components of the flow. As shown in Ricot et al. 


2009 the large wavenumber dissipation is greatly enhanced in order to stabilize 


the numerical scheme. All the Im(wj) in the Von Neumann stability analysis 
are greatly decreased after the cutoff value which decreases the accuracy of the 
propagation of high wavenumber components of the flow. In the present model 
only the dissipation of uig is increased and uj± are left untouched with respect 
to the LBM-BGK scheme which should provide a better accuracy for acoustic 
applications. 


4 Benchmarks 

In order to validate the model we are going to study a 2D and a 3D case, 
namely the dipole-wall interaction, and the turbulent jet. Both these flows are 
challenging from the numerical point of view since they exhibit a turbulent 
behavior (2D as well as 3D turbulence). 


4.1 Dipole—wall interaction 


This benchmark is based on the works of Clercx and Bruneau 2006 and Latt 


and Chopard 2007 . It analyzes the time evolution of a self-propelled dipole 


confined within a 2D box. The geometry of the box is a square domain [—A, A] x 
[—A, A], surrounded by no-slip walls. The initial condition describes two counter¬ 
rotating monopoles, one with positive core vorticity at the position {xi,yi) and 
the other one with negative core vorticity at (x 2 ,y 2 )- This is obtained with 
an initial velocity field Uq = {ux,Uy) which reads as follows in dimensionless 
variables 


U. = llwell ( 2 / - ^ llwell {y - (72) 

'^y = +\ ll‘^e|l ^ ||a;e|| (a; - 0 : 2 ) 6 “^’'''/’'°^". (73) 


Here, Xi = \/{x — XiY + {y — yiY, defines the distance to the monopole centers. 
The parameter tq labels the diameter of a monopole and We its core vorticity. 

The quantity we are interested in monitoring is the average enstrophy which 
is defined by 

= J oj'^{x,t)dxdy, (74) 


where to = d^Uy — dyUx is the flow vorticity. 

Under the actions of viscous forces, the dipole described by Eqs. (721 and (73) 
develops a net momentum in the positive x-direction and is self-propelled to¬ 
wards the right wall. The collision between the dipole and the wall is character¬ 
ized by a 2D turbulent dynamics where the wall acts as a source of small-scale 
vortices that originate from detached boundary layers. After the first collision 
the monopoles under the action of viscosity are re-propelled against the wall. 
These collisions give rise to two peaks of enstrophy (see Fig. |^. The value 
of these local maxima will be used for comparison with the results obtained 

2006] . Several snapshot of the 


with a spectral method in Clercx and Bruneau 


dynamics of the dipole-wall collision can be found on Fig. 
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Figure 6: Average energy (plain line) and average enstrophy (dotted line) evo¬ 
lution with time. The two enstrophy peaks (ili) and (^ 2 ) are highlighted. 



Figure 7: Snapshot of the vorticity at t = 0,0.15,0.0.32,0.4,0.48,0.64,0.72.0.8 
from left to right and top to bottom. Black is for positive and white for negative 
vorticity. 
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Figure 8: Numerical accuracy in the 2D dipole-wall collision flow for the two 
enstrophy peaks, (left) wi = 3313, and (right) uj 2 = 1418. The error curves for 
the enstrophy peak obtained with the present scheme and the MRT scheme. 


In this benchmark the initial core vorticity of the monopoles is fixed to 
We = 299.5286. Furthermore, the Reynolds number and the monopole radius 
are set to Re = LUjv = 2500 and rg = 0.1. The positions of the monopole 


centers are {xi,yi) = (0,0.1) and ( 12 , 1 / 2 ) = (0, —0.1). The approach of |Latt 
and Chopard |2007| is used to set up the initial condition. 


The error on the value of the enstrophy peak is the principal quantity of 
interest here. It is defined as 


f;, = |(fi,)-(r!,,ib)|/(ri*), 


(75) 


where f = 1, 2 is the label of the enstrophy peak. The value of the enstrophy 
computed with the LBM (either the present model or the MRT model) is noted 
The reference value (17^) is the value found in Clercx and Bruneau 
[2006| and is given by (Di) = 3313 and (172) = 1418. The convergence study is 
performed by keeping u constant and modifying I/ib while varying the resolution 
N (Re = UihN/u). Here C/ib = O.OliV/125 with N = 125,250,500,1000. This 
rescaling of the velocity has as an effect to reduce the compressibility errors 
from the simulation (see Latt ( 2 ^). As shown by Figs. the difference in 
accuracy between methods is not dramatically different. The differences appear 
when one pushes the numerical scheme to more challenging Reynolds and Mach 
numbers. 

In order to test the numerical stability and the ability to go to “higher” Mach 
numbers (but still lower than one) we also simulated the dipole at a maximal 
Ma number of 0.7 (corresponding to a characteristic velocity of C/ib = 0.032) 
and Re = 2500. At such high Mach number the MRT model was numerically 
unstable. The maximal stable reachable Mach number was of 0.46 (correspond¬ 
ing to characteristic velocity of U\h = 0.02). For this test the Mach number is 
kept constant and therefore one modifies the viscosity (in order to keep Re con¬ 
stant). By increasing the the resolution we do not remove the compressibility 
error terms (as discussed in Latt 2007 ). This explains the first order decrease 


of the error observed in Fig. |8[and the lower accuracy of the results. 

We notice that in this case the accuracy is much lower since the compress¬ 
ibility effects are much higher. Nevertheless the stability of the present model 
is highly enhanced with respect to the MRT model. Since such a “high” Mach 
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number was not achievable with the MRT model. 


4.2 Turbulent jet 

In this section we will perform the simulation of a turbulent round jet (see|Pop'e 


2005 ) at Re = NU/v = 6000 with N and U being the diameter and the speed 
of the jet respectively. The aim will be to recover correctly the self similar 
behavior and the correct energy spectrum, and also correct pressure spectrum. 
The computational domain is depicted on Fig. The domain size in units of 
the diameter of the jet was chosen to be of [50,30,30] x N. In order to avoid 


Sagau 

2013 

et al. 

2003 f 


Sagaut 2013 were added in the domain. Furthermore a vortex ring (see Bogey 


for example) is added at one diameter from the inlet to help the 
onset of the instability and allow for the development of turbulence in the flow. 
The Mach number of the flow is set to 0.4. The value is chosen to be rather 
large in order to really challenge the numerical accuracy and stability of the 
models. 

The quantities of interest are defined as follows. The velocity, u{x,r,9,t), 
is given in cylindrical coordinates centered around the center of the jet. The 
mean axial velocity field at the center of the flow is given by 


Uc{x) = {u^{x,0,0,t)), (76) 

where (•) is the time average. The jet’s half width, ri/ 2 {x) is defined such that 


{ua;{x,ri/ 2 ix), 0 ,t)) = -Ucix). 


We will also study the Reynolds stresses where 


U = U — (U) 


(77) 


(78) 


The simulation is performed with a very low resolution of = 10 points in 
the diameter of the jet and with no explicit turbulence model for the case of 
the present model. For the MRT case a Smagorinsky model was needed (see 


Krafczyk et al. 2003 , Malaspinas and Sagaut 2012 ) in order to obtain stable 


results. The fact that no explicit turbulence model is needed for our novel 


20 








































Figure 10: The jet’s half width with respect to the position for the MRT and 
the present model. The spread rate is of respectively of SruKT = 0.074 and 
S'rpresent = 0.093 for the MRT and present model. 



Figure 11: Non-dimensional velocity profile with respect to the rescaled position 
for the MRT (left) and the present (right) models. Five different x positions in 
part of the flow where the fluid is in a turbulent regime are depicted. 


scheme seems to imply that the regularization operation has the effect of an 
implicit turbulence model and would deserve a more in-depth analysis. 

depicts ri/ 2 ix) from which one can compute the spread rate Sr = 


Fig. 


10 


dri/ 2 /dx of the jet for the MRT and the present model. Although both models 
exhibit a self-similar behavior since there is a linear growth of the value of 
the spread rates are significantly different. One has respectively = 0.078 

and S'rpresent = 0.093. The expected value of the spread rate is of roughly 
0.1 (see Pope 2005|). Therefore the present model seems to provide a more 


accurate result than the MRT model. 

As shown in Fig. E] the self-similar behavior is observed for both models as 
for five different positions in the direction of the jet, the normalized average ve¬ 
locity profiles overlap for the MRT and for the present model. For the Reynolds 
stresses {{u'J'), and (u(,^) respectively) which are depicted on Figs. 


‘i 


one can notice that the self similar behavior is shown for the present mO' 

For the MRT model while close to the jet center the results seem self-similar 
(and also are coherent with what is observed with the present model), one can 
see that when going to r/ri /2 > 1.5 the Reynolds stresses are not overlapping 
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Figure 12: Non-dimensional Reynolds stress component with respect 

to the rescaled position for the MRT (left) and the present (right) models. Five 
different x positions in part of the flow where the fluid is in a turbulent regime 
are depicted. 



Figure 13: Non-dimensional Reynolds stress component (Wy^)/rtc with respect 
to the rescaled position for the MRT (left) and the present (right) models. Five 
different x positions in part of the flow where the fluid is in a turbulent regime 
are depicted. 


anymore and even worse, they are not converging towards zero as they should. 
This behavior can be explained by looking at the instantaneous velocity field 
at a given time. As one can see from Figs. 15 and (which represent respec¬ 
tively instantaneous snapshots of the velocity norm and of the pressure fields) 
the results obtained with the present model are far less noisy and no spurious 
modes can be observed. The only “spurious” modes present in the present reg¬ 
ularized model are due to the vortex ring used to trigger faster the transition 
to turbulence as seen on Fig. For the MRT spurious modes can be observed 


in the pressure field although this model is expected (see Lallemand and Luo 
Xu and Sagaut 2011] for example) to dissipate the pressure waves at a 


2000 


higher rate than for BGK models. 

Finally we also computed the energy and pressure power spectrum. One can 
see that for both the MRT and present model the —5/3 slope in the inertial range 
is recovered for the energy spectrum (see Fig. 17). For the pressure spectrum 


one expects to find a —7/3 slope in the inertial range as shown in George et al. 
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Figure 14: Non-dimensional Reynolds stress component with respect 

to the rescaled position for the MRT (left) and the present (right) models. Five 
different x positions in part of the flow where the fluid is in a turbulent regime 
are depicted. 
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Figure 15: Instantaneous velocity norm for the turbulent jet for the present 
model (left) and the MRT model (right) in logscale. The sponge zone region is 
removed from these pictures. 
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Figure 16: Pressure fluctuations for the turbulent jet for the present model (left) 
and the MRT model (right) in logscale. The colorbar scale has been reduced 
to allow seeing the small acoustic perturbations. The sponge zone region is 
removed from these pictures. 



Figure 17: Energy power spectrum with respect to for the MRT model (left) 
and the present model (right). 


1984 . While for the present model the pressure spectrum slope is correct, in 
the MRT case the slope of the pressure is closer to —5/3 (see Fig. 18). This 
difference indicates that our new model represents the dynamics of the flow with 
a much greater accuracy. 


5 Conclusion 

In this paper we demonstrated the existence of a recursive formula that al¬ 
lows the reconstruction of the non-equilibrium moments of the Boltzmann-BGK 
equation at any order, only by knowing the lower order moments. This prop¬ 
erty allows to propose a regularization procedure for the BGK lattice Boltzmann 
method that is increasing the overall accuracy of the method and removing the 
majority of the visible spurious modes present in standard BGK and MRT mod¬ 
els. This is shown by two benchmarks: the wall-dipole interaction (2D case) 
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Figure 18: Pressure power spectrum with respect to the frequency for the MRT 
model (left) and the present model (right). 


and the turbulent jet (3D case). Although in the 2D case the increase in ac¬ 
curacy of the model is not spectacular, the stability is greatly improved. For 
the 3D case a great improvement is shown. Not only the Reynolds stresses are 
found to be more accurately represented (with respect to the MRT model) but 
the pressure spectrum has the expected behavior. 

Finally one can notice that the model acts like an implicit large eddy sim¬ 
ulation model, since the turbulent behavior of the jet is reproduced accurately 
without any need to use an explicit turbulence model. Its relative simplicity and 
low cost of implementation make it very appealing for high Reynolds number 
and moderately high Mach numbers (of a maximal value of roughly 0.5). 

An interesting optimization of the present scheme could be achieved by using 
a D3Q19/D3Q15 quadrature instead of the D3Q27 used here. To do so the re¬ 
cursive formula used throughout this paper should be generalized to alternative 
basis vectors that are not Hermite polynomials. 

Finally the present approach may also be extensible for higher order quadra¬ 
tures (and higher order physics where one includes temperature for example) 
like thermal and compressible flows where one could reduce dramatically the 
memory needs of such cases (which is excessive for the moment since one needs 
D3Q121 quadratures). It could also provide a way to deal with boundary condi¬ 
tions for these kind of models, since with a very limited amount of information, 
one can reconstruct the complete populations with an accuracy consistent with 
the model. 
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A Computation of the recursive relation for off- 
equilibrium Hermite coefficients 


The aim of this section is to prove Eq. (311 that we recall here 

(") ("-1) , f (2) , ^ 

+ perm(a„)J . 

To prove it we need three relations. The first is the recursive relation of the 


equilibrium distribution Hermite coefficients (see Shan et al. 2006 , Malaspinas 


2009 ) 


7/ and - n 


(79) 


The second is that the order zero Chapman-Enskog expansion of the continuous 
Boltzmann-BGK equation leads to the Euler equations 


dtp + V ■ (pit) = 0, 

pdtu + pu ■ Vu = —Vp. 


(80) 

(81) 


And finally the order one Chapman-Enskog expansion of the Hermite coeffi¬ 
cients 


1 




— /f) -L r) 


(n-l-1) 


(82) 


where “perm” stands for all the cyclic permutation of indexes ai, 

Replacing n by n — 1 in this last equation and multiplying the result by 
one gets 

--91 —91 4 - 7 / ,9 

+ Ua„ (5aiao"a2^.La„_2 + pemi) • (83) 

Using the chain rule one can rewrite this equation as 

^ <^1 r», n, W “O n, /v. , , \ Uq,, I 


= 5tao”2,+'9a„+iao?J'i!Lc-„+i + ('9aiao"a2!.^..,a„ +perm) 


(i) 

' -V-' 

(ii) 

- (4”a2^.La„-2^“i“«/. +Pe™(a„)), 

^-V-^ 

{iii) 

(84) 

where perm(Q;„) is the cyclic permutation of all the indexes not equal to a„ (the 
a„ index never changes position). The (i) part of the above equation is equal 
to (see Eq. 


(i) = 


(85) 
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The (ii) can be rewritten 


(a) = {dtUa^ +Ma„+i9a„_^i?iQ„) 

- ^ai ' ' ' '^Ctn — i^anP^ 


( 86 ) 


where in the first equation we used Eq. (79) and in the second we used Eq. (811. 
Finally the {in) part reads 

{in) = -2 + perm(a„)^ - Ua-, ■ ■ ■ Ua^_^da^p, 

= ^ (ao?c7i^La.-2«pL-ia„ + perm(a„)) - ■ ■ ■ Ua„_,da^p, (87) 


where we used the chain rule in the first equation and Eq. (25) in the second 
equation. 

Finally adding {i), {ii), and {Hi) one obtains 


An) 


Jn-l) 




(uai-.. 


( 2 ) 


'^an-2^1,arv-ia„ 


perm(a„ 


B The higher order off-equilibrium Hermite co¬ 
efficients and implementation formulas 


In two dimensions the off-equilibrium Hermite coefficients a)" are given by 

l,yy -r I fJUyWyUx, ( 88 ) 


a^llyy = Uxa^lL + + TpuldyUx 


A3) 


I%xy = + Uya'flx + TpuldxUy, 


^l,xxyy ^y^xxy ^x^l,yy + 2uxUya[^l,y + rpUxUy {uldxUy + u^dynx) ■ (90) 


,(2) 


(89) 


In three dimensions the off-equilibrium coefficients at order n = 3 read 

O-llxy = ‘2‘Uxafly + UyOil,^ A TpuldxUy, 

* 

(3) o (2) I (2) , 3 O 

O-l'xxz = ^'U-xU)’ + U:,a\’ + TpUxdxU,, 


a?lyy = UxU^'i^yy + 2uya‘fly + rpuldyUx, 


o-iLz = Uxa\ i„ + ^ Tpuld.Ux 


A3) 


ilyzz = + rpuld^Uy, 

'' - ^ - ' 

* 

(3) o (2) , (2) , 3 a 

0\ yyz = ^Uya\ A- u^a\ + Tpu dyU 


( 2 ) 


y'-'y'Xz^ 


(3) _ (2) , (2) , (2) 

^l,xyz ~ '^x^i yz + UyUi ^j.^ + UzUi ^^y, 


(91) 

(92) 

(93) 

(94) 

(95) 

(96) 

(97) 


27 











while the n = 4 are given by 


^l,xxyy ^x^l,yy 

+ prUxUy 1 

(2) 1 (3) 

xUyCll xy 4" Uytti xxy 

{^UydyVjx 4“ ^x^X^y^ 

o-^ilxzz = + 2m 

4- prUxUz 1 

(2) , (3) 

x'^zO'l^xz 4“ '^zCll^xxz 

[^UxdxUz 4“ 2u^dzUx') 

„(4) _ 2 (2) 2u 

^l,yyzz — ^y^l,zz ^ 

+ prUyUz 1 

* 

(2) 1 (3) 

yUzal^^z: + Uza\^^yz 

^^y^y'^z 4“ 2uz;0zUy'j , 


* 


a 


(4) 

l,xyzz 


a 


(4) 

l,xyyz 


(2) , (2) , (2) , (3) 

UxUya\'^^ + + UyU^a\!^^ + 

+ prul {uxdzUy + UydzUx), 

'■ -V-' 

* 

O (2) , 2 (2) , (3) 

2u^Uya\!y^ + 

+ pTUxU^dyUz, 

' V-^ 

* 


,(4) 


,2^(2) 

'i.y. 

pTUyU"OxU^ 


^'l,xxyz ~ '^x^l,yz + 

.,3q 

r>^X 


(2) I (3) 


(98) 

(99) 

( 100 ) 

( 101 ) 

( 102 ) 

(103) 


(5) 
n ' 
l,xxyzz 


2 (2) I 2 (2) 

y dz^' 


'^xUya'^;^^ + U^U:,aiy^ + ZU^UyU:,ai, 
+ prUxUz {uluydxUz + 2uyuldzUx 

'' -V- 


M^UyU^ai;^^ 

-\- UxU^dzUy^ 1 


* 

~2v?u -^2u -^u 

^l.xxyyz ~ ^^x^y^l,yz ' ^^x^y^i^xz ' ^z^l,xxyy 
-\- pTU^Uy {Ux^x^z “1“ UydyUz\ 


sit 

(2) I o (2) I 2 (2) I (4) 

alL + ^UxUyUza\!y^ + UyUza\’xz + Uza\^'xyyz 

^'^Uyliz ip^y^z:^x 4 ” ^x^y^y^z 4 ” ‘2yUx'^z^z'^y^ i 
■ -/ 


(5) 2 {y^J 

^l,xyyzz ~ '^xUyO-i zz - „ - „ 

4“ P'^UyUz (^UyU^dz^X 4“ Ux^ydy 

'• -V- 


(6) 

n ' 

l,xxyyzz 


(104) 


(105) 


(106) 


* 

+ ‘^uluyUza^ilz + ‘^Uxuluza{]xz + ^zaflxyyz 
4“ P'^UxUyUz (^UxUydxUz 4“ 2UyU^dzUx 4“ UxUy^yUz 4“ 2UxU^dz^y'^ • 

(107) 


The * terms are the ones that are not present in the continuous case (see 
Eq. @) and are due to the discretization of the microscopic velocity space. One 


can see that as pointed out in Subsec. 3.3 these terms are of order 0(Ma"^ ) for 


the coefficients of order n, which makes them one order of magnitude smaller 
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than the other composing them. Therefore they are simply ignored for the 
computation of the terms. 
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